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^ ; Abstract 

For years, astrophysicists, plasma fusion and fluid physicists have puzzled over Rayleigh- Taylor 
. turbulent mixing layers. In particular, strong discrepancies in the growth rates have been observed 

between experiments and numerical simulations. Although two phenomenological mechanisms 
(mode-coupling and mode-competition) have brought some insight on these differences, convincing 

V ' 
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! analytical expression of the growth rate derived from first principles. It is compatible with both 

O ' mechanisms and is valide for a self-similar, low Atwood Rayleigh- Taylor turbulent mixing subjected 

C/3 ■ 

J>^' to a constant or time- varying acceleration. The key step in this work is the introduction of foliated 

Qh< averages and foliated turbulent spectra highlighted in our three dimensional numerical simulations. 

, We show that the exact value of the Rayleigh- Taylor growth rate not only depends upon the 

' acceleration history but is also bound to the power-law exponent of the foliated spectra at large 

vn ■ 

lO ■ scales. 

o 

^' 
o 

o\ 
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theoretical arguments are missing to explain the observed values. In this paper, we provide an 



PACS numbers: 52.35.Py, 47.20.Bp, 47.27.eb, 47.27.te 
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I. INTRODUCTION 



If allowed to evolve long enough, the Reynolds number of a Rayleigh- Taylor (RT) flow 
2| increases and the mixing eventually becomes fully turbulent 3|, Then, at late time, 
the width of the mixing layer grows according to 



L{t) = apAg{t)t'^ 



(1) 

when the gravitational field history [5| is of the form g{t) oc t^. The Atwood number A 
quantifies the density contrast of the initially pure fluids and is tied to the acceleration in 
any buoyant flow. From a purely theoretical point of view, the RT growth rate coefficient 
Op should be universal if the late stage of the evolution of the mixing zone is self similar and 
should only depend upon the exponent p of the gravitational field history. However, even 
today, the value of "alpha" (a), which is the name given by the "alpha-group collaboration" 
to ap/2 when p = (constant acceleration), is still a subject of controversy. 
For years, the discrepancy between its value inferred from numerical simulations 4, 6, 7| 
initialised with small wave lengths), scattered around 0.05, and its experimental values 
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111 ], twice as large, has thrown doubt on its universality. Not long ago, the heuristic 



mode-competition mechanism 
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13| put an end to this paradox: adding a small amount 



of modes of large wave length to the perturbation of the initial interface helps increase the 
value of a and reconciles numerical simulations with laboratory experiments. But, it is 
also instrumental in demonstrating that the initial conditions in any laboratory experiments 
do feature these modes of large wave length triggered by external causes like vibrations 



and residual motions. Flows of interest (geophysics, atmospheric, astrophysics [14, llSj), to 
a large extent, are not confined to a laboratory. Orders of magnitude may separate the 
length scales of the interface perturbation with the geometric size of the flow. This is why, 
the mode-coupling mechanism, initiated with small wave length modes, is thought to be 
realised in astrophysical mixing ^\ and other natural flows. 

Here, in order to find the theoretical value of the growth rate (op), we have adopted a 
point of view which has been used in the study of homogeneous turbulence (HT) for almost 
half a century but never investigated in RT turbulent mixing. The "large-scale structure of 



homogeneous turbulence" [16|,ll7| has been studied to understand the conjectured persistence 



of big eddies and the decay rate of various physical quantities in turbulent flows. At large 
scales, the turbulent kinetic energy density spectrum behaves as E{k) oc where the value 



of s depends on t 
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le initialisation of freely decaying HT in a way which is still unclear 
Bulk quantities of such flows, like turbulent kinetic energy, var y a s 



power-laws at late time and their exponents appear solely related to the value of s [23 1. 
That underlines the fact that large scale spectra have a strong influence on the behaviour 
of the overall flow. The theoretical study of freely decaying HT is made difficult by the 
fact that large scales are driven by non linear effect (triadic interaction). In contrast, the 
RT turbulence, although anisotropic, displays a significant difference with freely decaying 
HT. When plunged in a gravitational field, density fiuctuations produce motion through 



buoyancy. This production mechanism affects all scales [2J], even the largest, and produces 
turbulent kinetic energy at all times. It overcomes non linear transfer mechanism at large 
scales making the theoretical study of RT turbulence much easier. 

Here we demonstrate that it is possible to display an exact formula for ap which depends 
upon p, the exponent of the acceleration history, and s, the power-law exponent at large 
scales of foliated spectra that rise up naturally, along with foliated average, from the theo- 
retical developments hereafter. This foliation procedure is introduced for the first time in 
turbulent mixing theory because, among other interesting features, it allows to cancel out 
the formal effect of pressure in the equations, while keeping its physical effect on the whole 
mixing zone, making theoretical calculations amenable. 

II. TRANSVERSE AVERAGE 

A RT mixing occurs whenever two fluids of different density, initially separated by a sharp 
interface, are subjected to a gravitational fleld directed from the heavier (density ph) to the 
lighter (density pi) and when the initial interface is slightly distorted by a perturbation 
whose spectrum contains modes with sufficiently large wave lengths to overcome viscous 
effects at the onset. Such a flow, with two incompressible fluids in the low Atwood limit 



A = {ph — Pi)/ {Ph + Pi) ^ 1 (Boussinesq approximation), is governed |25| by a concentration 



equation ([2]), the Navier Stokes equation supplemented with a buoyant source term and 
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the incompressibility constraint (jlj) 

dtC+{uV)c = kAc, (2) 
dtu+{uV)u = -VP + 2Agc + u Au, (3) 
V -u = 0, (4) 

where c [x, t) is the mass fraction of the heavy fluid (c = in pure hght fluid and c = 1 in 
pure heavy fluid), u{x,t) the velocity field and p{x,t) the pressure field. The molecular 
diffusivity n and the kinematic viscosity of the mixture u are both assumed to be constant. 
The gravitational acceleration vector g{t) points down and its intensity may be time depen- 
dent as outlined earlier. In our simulations (described in the section §VIip . the flow fills a 
rectangular box of size H x H x 2H (see Fig{T]for more details). 

The calculation starts by averaging ([2l|4]). The transverse spatial average of a physical 
quantity q at height z is defined according to 

q{z,t) = ^ j ^ dxdyq{x,y,z,t) , (5) 

along with its fluctuating part q' {x, y, z,t) = q [x, y, z,t) — q [z, t). Velocity, concentration 
and pressure fields can be decomposed into a mean and fluctuating part: u = u + u', 
c = c + c' and p = p + p'. The velocity along the z axis will be called Uz and v will denote 
the transverse component in the xy plane so that u = {v,Uz). The relation ul{z,t) = is 
inferred from the divergence equation (jlj) and the fact that vertical velocity must be zero on 
the top and bottom walls. Homogeneity and isotropy in the transverse plane (no privileged 



direction) give = 0, v[c' = and v'^u'^ = at all times (25[ {i G {x,y}). The behaviour 
of the fluctuating flow will be studied at large scales (low k) and at high Reynolds number. 
Viscous and diffusive terms (in factor of u and k) can then be discarded. Therefore, the 
resulting set of equations governing the evolution of the fluctuating flow at large scales reads 

dtc' + u'.dzc + di {v[ c) + dz {u^ c) = 



dzu'^ c', (6) 
d,u'^ + d, {v'^ «0 + dz (K <) - dz iKK) = 

-dzp'-2Agc', (7) 



dtv: + d,{v'^v';} + dziu'A)-d,{v'^v';} = 

-d,p', (8) 



where Einstein summation convention is used on the repeated indices i and j G {x,y}. 



III. FOLIATED AVERAGE 

The last equations can be Fourier transformed for the purpose of showing that the growth 
rate of a RT mixing zone depends upon the structure of turbulence at large scales. Since the 
physical domain is assumed to be periodic along x and y (see Fig{T]), one defines a Fourier 
transform operator acting on these two directions only (transverse directions). It turns a 
generic quantity / (r, z, t), where r = {x, y) is the transverse position, into / (fc, z, t), where 
k = {kx, ky) is the transverse wave vector. Equations ([MH]) become 



dtd = —u'^ dzC — dziu'^* d) 

-tkjv^*?, (9) 
dfu'^ = —dz{u'^* u'^) — dzp' — 2Agc' 

—ikjv'-*u'^, (10) 
dtv[ = -dz{u'^* v-) - ikip' 

-ikjv'j*v'i, (11) 



where * is the usual folded product. In order to simplify flQl lTTl) one can eliminate the effect of 
dz, and with it, the effect of pressure and non locality, by integrating these equations along 
the anisotropic direction z. Therefore, these arguments suggest introducing the foliated 
spatial average which, for a physical quantity q, is defined according to 

{q)ix,y,t) = j Qi^^y^z,t) dz . (12) 

It amounts to chopping the domain (foliation) in elementary slices, transverse to the 
anisotropic direction, and to squeezing them by adding them up, resulting in an effective 
2D flow (see Fig|2]). Both in the concentration equation iQ and in the equation of motion 
( ITO|) . {dziu'^* q')) =0 (where q' = c' or u'^) because at top and bottom walls u'^ = 0. The 
average of the pressure gradient is {d^p') = because pressure fluctuations decrease outside 
the mixing zone. It has to be assumed that the top and bottom walls are far enough so 
that lim^^-|-ooP'(^; k,t) = 0. Consequently, the boundaries do not affect the flow within the 
mixing zone if L{t) <^2H (which can be relaxed to L{t) ^ H in simulations). This is why. 
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the resulting concentration equation and equation of motion after foliated average are 

dt{c') = -{u',d,c)-ik,{v'^*?) , (13) 
dt {uQ = -2A g (?) - % kj {v'j * M^) . (14) 

In eq. (fT3l) . the term dj:) can be simphfied because the profile of c{z, t) is a straight line 
(constant slope), from to 1, from one edge of the mixing zone to the other, at small Atwood 



ll| and numerically, in our simulations 



number. This is corroborated experimentally [lO, 
{A = 0.1), to within statistical fluctuations due, at a given height z, to the finite number 
N(t) ^ H'^/i{t)'^ of eddies in the transverse plane of area H^. The integral length scale i{t), 
i.e. the diameter of a typical eddy, varies like L(t) oc t^"*"^ in the self-similar regime. These 
statistical fluctuations become negligible in the limit H +oo since their rms amplitude is 
of the order 1/ ^yN{t) ^ i{t)/H. In this limit, the derivative dzC is then uniform within the 
mixing zone. Thus, the theoretical development suggests defining L such that d^c = 1/L{t) 
which can be brought close to another definition of the width at low Atwood number (see 
the end of §VIII) . Therefore, it is possible to write (u'^dzc) ~ {u'^) dzC = {uQ/L{t). As 
a consequence, (i) at large scale {k ^ 27r/£(t)), (ii) at low Atwood (iii) at high Reynolds 
number and (iv) in the limit H +oo (i.e. L{t) <^ H), the evolution of foliated second 
moments can be derived. From eqs. fll3tll4|) . straightforward algebraic manipulations provide 

(15) 

~ ~ ~ ~ (u')(u')* 

9,(Re[«)(cr]) = -2Ag{c'){cr-^-^l^ 

-k,^\ (17) 

where C| = 2Re[i{v'j * c')(c')*], ^| = 2Re[i{v'^ * <)(<)*] and = Re[z(t^. * <)(?)*] + 
Im[?(t'^- * c'){u'.^*] cancelled out when ensemble averaged. 

IV. ENSEMBLE AVERAGED FOLIATED SPECTRA 

Let us note q the ensemble average of a generic quantity q. Because of isotropy in the 
transverse direction in the self-similar fully turbulent regime, ^j{k,t) = 0, ^j{k,t) = and 
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^j^{k,t) = 0. It suggests defining the ensemble averaged foliated spectra for kinetic energy, 
concentration and production respectively as 

S,{k, t) = [ dQ , (18) 



2 „ 

£c{k, = k J dn {?){?*) , (19) 

suk, t) = j d^M{K){c*)] , (20) 

where the integration over the solid angle dQ is to be understood as the integration over all 
directions of k (in the transverse plane). It is then found that these spectra, referred to as 
foliated spectra for short, are coupled to each other at large scales according to 

4^2||, (21) 

4 = 2AgS,, , (22) 
= 2Ag£c + 2-^ . (23) 

L[t) 

It should be stressed that this closed set of equations does not hold if foliated average 
operators are removed or even replaced by transverse average because of the non local 
effect of pressure. Since these equations are linear, autonomous in k and do not display 
derivatives with respect to k, it comes as a byproduct prediction of the theory that, with 
the aforementioned definitions, all three foliated spectra must have the same power-law 
exponent at large scales: k^. The exponent s is independent of time to comply with self- 
similarity. This common value is a special feature of the foliated spectra (see Fig|3D. Foliated 
spectra have deliberately been designated by script S in order to stress the difference with 
transverse spectra classically defined by 

E,{k,z,t) = ly dnu'.uC , (24) 

Ecik, z,t) = k j dVtdc* , (25) 

Ec.(fc, z, t) = -k j dnRe[<?*] . (26) 

Both types of spectrum were compared to show that this common exponent is specific to the 
foliated spectra whereas transverse spectra (concentration and velocity for instance) have 



different power-law exponents j?! at large scales. This can be understood when realising 
that a transverse spectrum is calculated using Fourier modes on a slice at a height z in the 



middle of the mixing zone. A foliated spectrum, on the other hand, uses the sum along z of 
all these modes. The resulting spectrum is radically different and benefits from interferences 
between modes at different heights. 



V. GROWTH RATE FORMULA 

For self-similarity to hold, the late time behaviour of the foliated spectra must be 

S,{k,t) = S'^kH'^ , (27) 
S,{k,t) = S^k't^'^ , (28) 
£,,{k,t) = E^kH^- , (29) 

at low k where £^ and £^.^ are constants independent from k and t. In addition, the 
exponent of k in each of this averaged spectrum (s) and the time exponents (cc, e^,, Ccz) 
must be constant in time. When replacing these expressions with Ag{t) = Agpt^ and 
L(t) = apAgpt^^"^ in eqs. ( !2Tll23|) one gets the following relations 

= e, + 2{p + 1) , (30) 

ecz = Cc+ip + l) , (31) 

£° = {ap{AgpfeJe,)S',, (32) 

£i = {apAgpej2)8l , (33) 

oip = , (34) 

^z 

which are not modelled but come straight from first principles. Using Eqs.( !27ll29|) and 
eqs. (130113^ ■ a non trivial exact result can be deduced: Sczi^^ ^)/ {'^^z{k, t)Sc{k, t)) = 1 
valid at large scale (see FigJl]). 

A fundamental feature of concentration is that its fluctuation variance, c'c' tends to a 
constant in the self-similar regime because concentration is a physical quantity bounded 
between and 1. It cannot go to zero for that would mean the mixing tends to be heteroge- 
neous in the middle of the mixing zone, thereby going against the fact that at all times any 
side of the mixing zone must be supplied with pure fluid from the other side. The mixed- 
fiuid-supplying-channels have been highlighted in FigH] and Figj2l This is exemplified by 



numerica 
to 0.8 



simu. 



0, Q, Q, 



ations and experiments reporting molecular mixing rates ranging from 0.7 
25|. 



The question now is: how about {c'){c')7 A very simple argument can be given: 
the discretized fohated average (c') = dz/H ^fjt]{/dz '^^^ ^le approximated by (c') ^ 
dz/H X]i^*-L(t)/d2 '^i since c' is only non zero within the mixing zone. That is why, {c'){c') ~ 

2^i=-L{t)/dz 2^j=-Lit)/dz ^i^i ^^'^ ~ jU 2^i=_mydz 2^j=-L{t)/dz ^0 ^ U " 

we say c'c' = Cg in the self similar regime as described in the previous paragraph. Therefore, 



(c')(c') ^ dz^H^ Y.\=ltlt)idz Co = dz/H'^ L{t)cl. As a result (c')(c') must grow like L{t) 
(corroborated by our numerical simulations). 



Therefore, since {c'){c') = J^^ dk Sc{k,t) varies as Jq^ dkEdk.t) oc 

^ec-(p+2)(s+i) must evolve as L{t) oc t^^"^^ the value of Cc is bound to be (p + 2)(s + 2) and, 
by inference, the time exponents = (p+2)(s + 2) + 2(p+l) and Ccz = (p+2)(s + 2) + (p+l). 
Therefore, the value of the Rayleigh- Taylor growth rate for the turbulent mixing zone width, 
ttp, is shown to depend upon which is the acceleration exponent, and s, the foliated spec- 
trum power-law exponent at large scales, according to 

g 

""^'^ ^ (p + 2)(s + 2)(p(. + 4) + 2(. + 3)) • ^^^^ 

This expression does not depend on Qp as expected and confirmed in our simulations. 

The demonstration leading to this formula did not provide any arguments against the 
idea of an s varying with p. However, growth rates, at various values of p given by our 



simulations in conjunction with those provided in [26j, were fitted to the theoretical formula 
and a remarkable collapse of the data was found for s = 4.0±0.1 assuming s was independent 
of p (see FiglS]). It is possible to check this result by using an other obvious method: direct 
inspection of the foliated spectra. The power-law at large scales of foliated spectra was 
checked on every simulation (p = 0, 1, 2 and 3) and was found to be compatible with 
s = 4. The formula (l35ll is therefore a predictive formula: it gives the value of ap knowing 
p (a controled parameter) and s or, the other way around, it gives s knowing p and ap. 
Our theory did not provide a value for s for it certainly depends upon the way the mixing 
flow reaches the self-similar regime. In order to find its value, the theory would have to 
cope with the transition regime to make the connexion with the linear growth of the initial 
conditions (s = 4 seems to be compatible with annular spectrum initial conditions). It 
must be emphasized that the exact result ( |35l) is a proof that, in the self-similar regime, 
the growth of the mixing zone only depends upon the structure of turbulence at large scales 
(appart from p, it depends on s only). 
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VI. IMPLICATIONS 



Obviously, the formula for ap carries information on the dynamics of the mixing zone 
width. Not so obvious is the fact that it also provides information on the way this dynamic 
affects the mixing at large scale. If it is assumed that the late time evolution of the mixing 
zone width does not depend upon initial conditions - which may be true when the charac- 
teristic length scale of the initial perturbations is much smaller than that of the physical 
domain if - it is reasonable to admit that it must depend on L(t) (and its derivatives) and 
g{t) which are the only control parameters left in the problem. From this assumption, basis 
of the buoyancy-drag approach j^, Q, it is possible to build a simple evolution equation, 

L{t)=C,Ag{t)L{t)-CdL^{t)/m , (36) 

which depends on two adjustable constants. The constant Cb quantifies buoyancy related to 
the mixing since Cb A can be seen as an effective Atwood number: the fluids are not pure, 
with density pi and ph, in the mixing zone but they are partly mixed with intermediate 
densities and the smaller Cb, the smaller the effective Atwood number, the stronger the 
mixing. The constant Cd quantifies drag, i.e. the exchange of momentum between raising 
and falling mixed-fluids structures. This phenomenological approach does not by itself 
provide any relation between these two adjustable constants. However, replacing the self- 
similar values ([I]) of L and g in ( l36l) allows to recover the exact formula ( l35l) if and only 
if Cb = 4:/{s + 2) and Cd = (s + 2)/2, which yields Cb = ^/Cd- This correspondence 
between the simple model and the exact result demonstrate the little importance of small 
scale as opposed to large scale in the bulk dynamics of the mixing. Furthermore, it enables 
to understand the influence of s on "intuitive" physical mechanisms like buoyancy, through 
the effective Atwood number, and drag. Accordingly, the bigger s, the stronger the drag, 
the smaller the effective Atwood number and the better the mixing. 

This is in agreement with numerical results 3, [l3| suggesting that a smaller proportion 
of long wave length in the spectrum (bigger s: mode-coupling) decreases the growth rate 
of the mixing zone width in opposition to a bigger proportion of long wave length (smaller 
s: mode-competition). Therefore, in the framework of foliated spectra where s is defined, 
the two phenomenological mechanisms of mode-coupling and mode-competition in Rayleigh- 
Taylor turbulent mixing can be explained and brought together. 
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There is no reason in principle why techniques developed in this work could not be 
applied to other mixing flows with one anisotropic direction. For instance, in the case of 



Richtmyer-Meshkov (RM) turbulence 27LAg = and the mixing width L(t) varies like 



in eqs. fl2m^ . Persistence of big eddies [I6j is then an exact result for the foliated RM flow 
since the rhs of eq. (1221) vanishes. The same type of reasoning that led to eq. fl35p allows to flnd 
6 = 2/{s + A) (also an exact result when s is the power-law exponent of ensemble averaged 
foliated spectra, at large scale) which for 3 < s < 4 predicts a 1/4 = 0.250 < < 0.286 = 2/7, 
in agreement with experimental and numerical results. 

VII. SIMULATIONS 

A set of four simulations of Rayleigh- Taylor turbulent mixing flows were carried out with 



the incompressible code SURFER [28|] for different acceleration histories g{t) oc t^ with 
p = 0, p=l,p = 2 and p = 3. For our purpose, gravity was added to the original version. 

An incompressible code is valuable to investigate flows under variable acceleration because 
hydrostatic pressure balances instantaneously with the acceleration fleld. On the contrary, a 
compressible simulation would have to cope with the creation of acoustic waves that would 
travel at the speed of sound to balance pressure with time-varying acceleration. These 
acoustic waves would bounce on the domain walls, would go back and forth through the 
simulation domain and would be detrimental to the simulations. 

Each of those simulations was performed with the same grid resolution: 128 x 128 x 256 
(along X, y and z respectively). The large scale outcomes of a numerical simulation of a 
Rayleigh- Taylor turbulent flow do not vary signiflcantly, for a given type of initial condition, 
by increasing the resolution. This explains the choice of such a modest resolution allowing 
to carry out more simulations. 

SURFER is parallelised, 3D and evolves two immiscible fluids separated by an interface. 
To reconstruct and advance the fluid interfaces in time, SURFER uses an exactly volume 
conserving variant of the Volume of Fluid algorithm with a Piecewise Linear Interface Cal- 



culation method (VOF/PLIC) 29|. The density of each fluid ph (for heavy) and pi (for 

light fluid) are constant in time. The Navier-Stokes equation governing the evolution of the 

velocity fleld u is given by : 

du 

p-^ = -Vp +V ■ {r]S) + arSsn + pg , (37) 
11 



where p is the density, p is the pressure and g is the acceleration. The dynamic viscosity rj 
equals rjh = phi^ or rji = piv (where is a common kinematic viscosity) and S is the rate of 
strain tensor defined by Sij = diUj + dj Ui. The surface tension a depends on the particular 
two fiuids that will be heterogeneously mixed and the delta function, 5^, is concentrated on 
the surface of the interface, r is the mean curvature of this surface and n is the unit normal 
on the surface. The discretisation of (137|) is performed on a MAC-type staggered grid and 
the pressure is computed using an iterative multigrid Poisson solver. 

The simulations were carried out in nondimensional units. Initially, two vertically stacked 
fiuid layers of different density, such that ^ = 0.1, have been considered (see FigH]). Periodic 
boundary conditions were prescribed on the four vertical domain walls, whereas no fiux and 
no slip conditions were imposed on the two horizontal walls of the domain, at the top and 
at the bottom. Initially, the velocity field u (a;, y, z,t = 0) is perturbed around the interface 
using a sum of random small amplitude modes which comply with the incompressibility 
condition. The wave numbers selected verify 15 < n < 17 (n ^ 1) to get a late time 
self-similar evolution of the fiow (mode-coupling mechanism). 

The value of ap is affected by how the width of the mixing zone, L{t), is measured 
or computed. Experimentally and numerically, two methods are employed and they both 
use c{z,t). The threshold method prescribes -Lth(^) =| -299 — Zi \ where c{zi,t) = 0.01 and 
0(2:99, t) = 0.99. This method is affected by statistical fluctuations of c producing noisy Lth(t) 
and an even noisier derivative. The integral method prescribes Lint(t) = 6/ dzc{z,t){l — 
c{z,t)) which provides an exact value if the proflle of c is affine. Since c is almost affine 
everywhere in our simulation at low Atwood, the difference between Lint(t) and 1/dzC is 
negligibly small. In a sense, the reasoning that led to the growth rate formula dictated the 
deflnition of L, by the integral method, which has been used to calculate the value of ap in 
each of the four simulations. 
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FIG. 1: Image of a turbulent Rayleigh- Taylor mixing zone for a constant acceleration {p = 0) when 
L{t) = H from our numerical simulations. The axes are defined in such a way that the mixing 
zone grows along z and is spatially periodic along x and y. The initial interface between pure light 
fluid and pure heavy fluid is located at z = 0. Thus, along the z axis at time t, the flow goes from 
a turbulent mixing of two fluids when | z |< L[t)/2 to an intermittent border at | z |~ L[t)/2 and 
finally to a laminar pure fluid (light or heavy depending upon the direction) when | -z |> L{t)/2. 
The vertical velocity u'^ of the flow is plotted in false colour on the vertical sections. It puts to the 
fore large structures going up (red) and down (blue), made up of mixed fluids. The 3D bubble-like 
shapes on top (at z ~ H/2) correspond to isosurfaces of at 80% of its maximum positive value. 
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FIG. 2: Image of thefoliated concentration and of the foliated vertical velocity. The top slice 
displays {c'){x, y) in false colour at the instant when L{t) = H for p = 0. The bottom slice displays 
— {u'^){x,y) (the minus sign is introduced so that colours match between top and bottom). On 
the colour chart, the value 1 corresponds to the normalised maximum. The striking similarity in 
location and shape of large scales on both planes (peculiar to foliated average) reveals that heavy 
mixed fluid tends to concentrate where the flow goes down (red) and light mixed fluid where the 
flow goes up (blue). 
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FIG. 3: Normalised foliated spectra versus transverse spectra. The left figure represents the ratio of 
velocity and concentration spectra, £z{k) / {£c{k) g{t) L{t)) in solid lines and Ez{k) / {Ec{k) g{t) L{t)) 
in dashed lines. The right figure represents the ratio of production and concentration spectra, 
£cz{k)/ (^£c{k) \l g{t) L{t)^ in solid lines and Ecz{k)/ (^Ec{k) \/ g{t) L{t)^ in dashed lines. Displayed 
results come from the simulation at constant acceleration p = 0. The colours correspond to different 
times in the evolution of the mixing zone: ti (blue)< t2 (red)< (green) defined by L{ti) = 0.5 H, 
L{t2) = 0.75 H and ^(ts) = H. In the whole domain of wave numbers, ratios of foliated spectra 
are distinctively smoother than transverse spectra. In the domain where n < 6 (large scales), the 
ratio of foliated spectra is constant and does not vary with time whereas, in the same domain, the 
ratio of transverse spectra can vary up to one decade. 
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FIG. 4: Normalised correlation spectra. The evolution of lZ{k, t) = 2g fkt)s\k t) lines) and 

of R{k, t) = 2 E (kt) E (k t) (dashed lines) at constant acceleration p = is compared at different 
times, ti (blue) < t2 (red) < (green). The simulation, at p = (but at p = 1, 2 and 3 as 
well), shows that when n < 6, 7?. = 1 to within 6% and with the same characteristic smoothness as 
already depicted on Figl2] (vertical axis is linear). On the contrary, R varies significantly as time 
goes by on the same range of wave numbers. 



20 



2 



o 

II 

CO 

o 

<x> 

4J a. 



3 

fl a. 



1.5 



0.5 



°2 



I i 



-10 12 3 

acceleration exponent: p 



FIG. 5: Observed to predicted ratio of growth rate ap. Green squares correspond to results 
obtained with SURFER and performed by the authors and red circles to results performed with 
TURM0IL3D by D. Youngs and A. Llor. A least square best fit has been applied on the numerical 
results using the formula (j35p and it was found, assuming s does not depend on p, that s = 4.0ib0.1. 
Error bars correspond to the size of the biggest interval containing all three different measures of 



1) ar 



L{t)/{Ag{t)t^), (2) a. 



L(t)/((p + 2) Ag{t) t) and (3) the same method described 



in ^. Error bars were not provided in 
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